Ion Sizes and Finite-Size Corrections for lonic-Solvation Free Energies 
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. Free energies of ionic solvation calculated from computer simulations exhibit a strong system size 

' dependence. We perform a finite-size analysis based on a dielectric-continuum model with periodic 

^\ \ boundary conditions. That analysis results in an estimate of the Born ion size. Remarkably, the 

finite-size correction applies to systems with only eight water molecules hydrating a sodium ion and 
results in an estimate of the Born radius of sodium that agrees with the experimental value. 
', Calculation of ionic-hydration free energies from computer simulations require great care|-to avoid artifacts due 
' to the long-range electrostatic interactions. We have recently shown that Ewald summationEJa can give jresults for 
single-ion free energies that are essentially system-size independent for as few as about 16 water molecules .cl That has 

■ been achieved by including the self-interactions O.Sg^^Ew of the ion (^ew ~ — 2.837297/L with L the length of the 
cubic box). Ewald summation determines the electrostatic interactions using lattice sums for a periodically replicated 
simulation box. This provides a natural description of the electrostatics in the periodic space resulting from the 
periodic boundary conditions commonly used in computer simulations. Deviations from the approximate finite-size 
correctionB'u are expected for soWents with finite dielectric constant e < oo (i.e., in the infinite-dilution limit of an 
ion in a non-conducting solvent),™ and if the ion size is comparable to the dimensions of the simulation box. 

^ , The self-interaction of an ion is introduced in the Ewald summation through interactions with the compensating 

■ charge background implicit in the Ewald summation approach .and with the peiriodic images of the ion. A similar 
O ' correction has been developed for free energies of polar molecules.Q Figueirido et al.H used a point ion in a homogeneous 
c/j dielectric medium as a model to extend the system-size corrections to solvents with finite dielectric constants e, finding 

_0 _ that 0.5q2^Ew,(l - 1/e) should be a good approximation to the finite-size effects. Here, we further extend Figueirido 
CZJ , et al.'s studyd to ions of finite size. We show that the ion-size correction is important for small system sizes. We also 
show that the analysis of the system-size dependence of ionic free energies leads to the definition of an ion size in 
'"^ , dielectric media in excellent agreement with the Born radius.D 

We approximate the free energy of charging an ion as the sum of the explicit simulation contribution Fgun and a 
self-interaction O.Sq^^Ew, as in Ref. ^. To account for further effects of finite system size, we add the difference of 
charging an ion in an infinite and finite system, F^o — Fl , 
> ■ 

5^ ■ F,i„, + ig2^E„ + (i^oo - Fl) (1) 

o : ^ 

The finite-size correction Fac — F^ is evaluated for a simplified model of a point charge g in a spherical cavity with 
radius R inside a dielectric continuum with dielectric constant e, as schematically shown in Fig. n. The finite system 
is formed by an ion in a periodically replicated box of length L, where the box is charge neutral through addition of a 
homogeneous background charge with density ~q/L^. The infinite system is obtained by taking the limit L ^ oo and 
Q i: corresponds to the Born modelB If the model used to calculate F^o — F^ were exact, F^ would cancel Fg^^ + O.Sg^^Ew 
However, the dielectric model is only an approximation to the simulation system with explicit solvent, making the 
finite-size correction only approximate. 
^ The finite-size correction is the difference of charging the ion in an infinite {L — > oo) and finite system, 

Foc-Fl = -q^ [e(e, R,L^^)~ ^(e, R, L)] , (2) 



^(e, R, L) is the electrostatic potential minus the bare self-interaction at the position of a unit charge, 

^(e,i?,L) = limf0(r)-l/r] . (3) 
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We determine the electrostatic potential (/)(r) for the dielectric model described above by solving the corresponding 
Poisson equation, 

V [e(r)V(/)(r)] = -Att [5{y) - l/L^] . (4) 
We define potentials 0-(r) and 0+(r) inside and outside the sphere. For R < L/2, the boundary conditions are 
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0_(r) = 4>+{r) for |r| = i? , 

'^-W..^ for in 
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The last condition reflects the periodicity of the potential and applies analogously for y and z. 

To solve this electrostatic problem, we expand and 0+ into a complete set of functions that satisfy Poisson's 
equation Eq. (||) and then choose the expansion coefficients to satisfy the boundMy conditions Eq. (^). Such functions 
are the kubic-harmonic polynomials Kn introduced by von der Lage and Bethe.13 We define 
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For a finite set of functions A'2„ up to polynomial order In < N , the solution is only approximate. We choose the 
coefficients in the least-square sense, defining a functional, 

= sine del [0-(r) - 0+(r)]' 
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with spherical polar coordinates r, c/j and 0. The last integral extends only over a circular region because this results 
in an analytical solutiaiLjo£_that integral. Also, the deviations from Eq. (pq) of truncated kubic-harmonic expansions 
of the Ewald potentialDa llj are largest in the corners of the box. Fitting only a circular region therefore produces 
smaller distortions of the potential at the center of the box. 

Minimizing with respect to the coefficients 6_i, 6oj 64|-etc. results in a set of linear equations for those 
coefficients. Using the computational algebra package REDUCEji3 we solve for the unknown expansion coefficients 
for different orders TV of the kubic-harmonic expansion. We then apply Taylor expansion with respect to 1 /e and the 
size parameter R/L. This results in an asymptotic expansion of i^(e, R, L) as 



e(e,i?,L) = 



Co 



eL 3L3 



^3g3 15^3 



(8) 



Notice that because of the fini 
in the simulation box is zero,!^^ 



size, 0(r) is normalized by adding a constant C such that the average potential 
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for 



00. The familiar Born term 0.5(1 — e ^)/R is thus contained in the normalization 



rather than 0(r 
constant C. 

The coefficients Cq, C2 and C4 are listed in Table |. As the order N of expansionpfunctions increases, the coefficient 
Co converges towards the self-energy of a point charge in a cubic Wigner lattice,ll3 cq L^Ew ~ —2.837297. The 
coefficients C2 and C4 appear to converge to values close to zero. The most interesting term in Eq. (^) is the lowest 
order correction for the finite size of an ion, —2n{e — l)i?^/(3eL^). Eq. (||) converges to the correct limit for e — > 00, 
which can be found independently using the continuity of 4>{r) at the dielectric interface |r| = R and the spherical 
geometry for e ^ 00, 
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(10) 



The last term is contained in the O (R^/L^, e"^) term of Eq. (^). 

We can now calculate the finite-size correction to the free energy of an ion with radius i? in a dielectric medium 
with Ewald-summation electrostatics, 

F^-FL = -^q^a^,R,L) , (11) 

using Eq. (||) and ^(e, R,L ^ oo) = for the particular choice of C [Eq. (^]. This results in an approximation to the 
free energy of charging the ion including finite-size effects: 

i^o « i^shn + Kew + ^ j . (12) 

In the limit e — > oo, the finite-size correction F^o — Fl does not contain terms of order and f^^- This explain^L 
the success of using Fgim + O.Sg^^Ew alone for the free energy of charging an ion in a conduct ingliil or highly polarti 
environment, without further finite-size correction that takes the ion size or the dielectric constant of the solvent into 
account. 

The validity of the approximate finite-size correction Eq. ( p^ is illustrated in Fig. |2|for a sodium ion in water. We 
fit the sodium-ion data Fsim for electrostatic free energies calculated from simulations with M = 8 to M = 256 water 
moleculesS to Eq. ( p^ ) with Fg and R as parameters. For the dielectric constant, e ^ oo is used, but values of e « 80 
have little effect on the result. That fit yields Fq = —402.5 ± 1 kJ mol~^ for the electrostatic contribution to the 
solvation free energy of sodium. The radius R of the sodium ion is found to be i? = 0.18 nm, injexcellent agreement 
with the effective Born radius 0.180-iim of that ion, as determined by Latimer, Pitzer and Slansky££l or, more recently, 
0.188 nm, as compiled by Marcus.llZl In the fit, the ion radius R is determined mostly by the data for small system 
sizes. Using only the M = 8 and 16 data and fixing Fq at —402.5 kJ mol~^, a fit of i? yields 0.159 nm {M = 16), 
0.190 nm (M = 8) and 0.184 nm (M = 8 and 16). 

The simulation result for the Born radius is 0.172 nm.0 That is, Fq and R are self-consistent assuming an uncertainty 
greater than 0.01 nm in R. We can therefore fit the free-energy data with the solvation free energy Fq as the only 
parameter by substituting the Born expression Fq « —0.5(7^(1 — e^^)/R for i?, 

Fq « + 2^ ^ (^&w + ) ■ (1^) 

That fit yields Fq ~ — 403± 1 kJ mol^^ (for e oo). These results show that the system-size effects on the hydration 
of sodium can be described accurately for M > 8 water molecules by the dielectric continuum model with periodic 
boundary conditions shown in Fig. ^ It is remarkable that continuum-electrostatics approximations apply even for 
as few as eight water molecules solvating the sodium ion provided the periodic boundary conditions are considered. 
Of course, determination of the value of Fq and the Born radius requires further molecular considerations. To-treat 
solutes with more complicated shapes will require additional computational effort, as was suggested previously.ES 

In summary, we have found an approximate finite-size correction for ions in water that takes into account the size 
of the,jjCtfis as well as the dicle&tpc constant of the solvent. These findings extend the previous work by Hummer 
et alffllB and Fi gueirido et aLBE^ The results of Refs. || and |ll| are recovered in the limit of a conducting solvent 
(e — > oo) and a point ion R/ L 0. In the limit of a point ion but for a finite dielectric constant e < oo, we obtain 
the finite-size correction of Ref. |^. For strongly polar solvents, such as water (e w 80), and typical system sizes of 
hundred ocjaiore solvent molecules with ions of comparable size, the difference to the finite-size correction obtained 
previouslyB'Q for e ^ oo and R/L — > will be small (< 1-2%). We have found that the finite-size correction derived 
in this paper gives qualitatively and quantitatively correct behavior. This finite-size analysis results in an estimate of 
the sodium-ion size that is in agreement with the experimental data. 
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TABLE L Coefficients of the approximation Eq. (|8|) to ^ for different levels A'^ of kubic-harmonic expansions. (A'^ = 2 means 
that no kubic harmonics have been used.) 



N 


Co 


C2 


C4 


2 


-2.9037 


0.9751 


-0.8581 


4 


-2.8309 


-0.1638 


0.1558 


6 


-2.8362 


-0.0728 


0.0711 


8 


-2.8398 


0.0371 


-0.0366 


10 


-2.8370 


0.0022 


-0.0023 



FIG. 1. Schematic representation of the contributions to the electrostatic solvation free energy of an ion, 

F.in, + O.bq^^E^ + (Foo - Fl). 

FIG. 2. Finite-size correction for the hydration free energy of a sodium ion. Shown is the free energy of charging a sodium 
ion in water from charge zero to e as a function of the inverse simulation box length 1/L. The symbols are simulation data 
for Ewald-summation electrostatics from Ref. ^ for M — S, 16, 32, 64, 128, and 256 water molecules. The dashed line is the 
fit to the correction formula of Ref. ^ (for M > 16), corresponding to Eq. ( p^ ) with e ^ oo and R = 0. That fit yields an 
extrapolated free energy Fq of —405 kj mol~^ . The solid line shows the fit to the correction formula Eq. (^2|) with e — » cxa, where 
the sodium-ion radius R and the free energy of charging Fo are estimated to be 0.18 nm and —402.5 kJ mol~^, respectively. 
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